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Abstract 

The Immersed Boundary (IB) method is a widely-used framework for the simulation 
^^ of fluid-structure interaction problems. The IB formulation maintains an Eulerian 

Q>^ , discretization for the fluid equations of motion while maintaining a Lagrangian 

r^ I representation of immersed elastic structures. The interaction between Lagrangian 

and Eulerian variables is mediated by convolution with delta function kernels. Our 
^^ . specific motivation is the modeling of platelets in hemodynamic flows. We explore 

(<— ^ I the use of a new geometric representation for the Lagrangian structures - radial 

cn ' basis functions (RBFs) - within the IB method. We compare our new RBF-IB 

method with the traditional IB method in terms of differences in Lagrangian marker 

positions, computational cost, maximum stable time-step size and volume loss. 

I 
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1 Introduction 



Intravascular blood clots (thrombi) are initiated by damage to the endothelial cell lining of a 
blood vessel and involve the formation on the damaged surface of clumps of cells intermixed 
with a fibrous protein gel. The cells involved in this process are platelets, and the subject of 
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this paper is a new approach to modehng platelets in order to simulate their adhesion to the 
injured vascular wall and cohesion with one another during the formation of a thrombus. As in 
our previous work, the Immersed Boundary (IB) method is used to describe the mechanical 
interactions among a collection of discrete platelets, the background fluid, and the vessel 
wall, but the representation of platelets within the IB method is very different than what 
has been used previously. Radial basis functions (RBFs) are now used in order to achieve 
more accurate and less costly simulations. 

In this introduction, we briefly describe the relevant biology, sketch the immersed boundary 
method which is central to this work, describe how the IB method has been used in our 
previous platelet aggregation simulations, and give an overview of how use of the new method 
changes this description. 



1.1 Biological background 



In the unactivated state in which platelets circulate, they are fairly rigid discoidal cells with 
a major axis of 2 to 3 microns and an aspect ratio of close to four. Unactivated platelets do 
not adhere to one another or to the intact endothelial cell lining of the blood vessels. While 
very numerous in the blood (f^ 250,000/yU.liter), they are much smaller and far less numerous 
than the red blood cells which comprise 40-45% of the blood's volume. It is the red blood 
cells rather than individual platelets that determine the blood's rheological properties [25]. 

Disruption of the endothelial cell lining exposes to the blood collagen and adsorbed von 
Willebrand factor (vWF) molecules in the sub endothelial matrix. Platelets adhere to both 
molecules via specific receptor molecules on the platelets' surfaces. In addition to slowing or 
stopping platelet motion over the subendothelium, this binding triggers intracellular signaling 
pathways that lead to platelet activation [17,23]. 

Among the many changes that an activating platelet undergoes, two are most relevant for 
this paper. One is that aubP^ integrin receptors embedded in its surface membrane become 
activated and capable of binding dimeric fibrinogen molecules and multimeric vWF molecules 
from the blood plasma. By binding to receptors on two platelets, these molecules serve as links 
between the platelets that allow them to cohere. The second is that the platelet's cytoskeleton 
is reorganized and as a result, the platelet initially becomes more spherical. However, it also 
becomes sufficiently flexible that over time it can spread out over the surface to which it 
is adhered. Activation can also be triggered when a platelet is exposed to sufficiently high 
concentrations of specific chemicals that are secreted into the plasma by other activated 
platelets. Thus, platelets near the injury (but which do not directly contact the damaged 
vascular wall) can be activated, and these chemically-activated platelets can also cohere with 
one another and with the wall-adherent platelets. Through these mechanisms, a platelet 
clump or aggregate grows at the injury site. As it grows, the aggregate can profoundly 
change the local fluid dynamics (to the extent that vessel occlusion can occur). Conversely, 
the changing fluid dynamics influence the further growth of the platelet aggregate both in 



terms of changes in the transport of platelets and chemicals to and from the injury site, and 
in terms of fluid forces which affect the binding and unbinding kinetics of the bonds between 
platelets and the subendothelium and between pairs of platelets. 



1.2 Modeling the mechanics of platelet aggregation 



Our platelet aggregation models [5, 8-10, 30] track the motion and behavior of a collection 
of individual platelets as they interact with the suspending fluid, one another, and the 
vessel walls. They also track fluid concentrations of platelet activating chemicals, cell-cell 
and cell-surface forces, fluid motion, and the local fluid forces on the growing thrombus. 
In the models, nonactivated platelets are activated by proximity to reactive sites on the 
injured wall, or through exposure to a sufficiently high concentration of activator in the 
fluid. Activation enables a platelet to cohere with other activated platelets, and to secrete 
additional activator. The platelets and the secreted chemical move by advection with the 
fluid and diffusion relative to it. Each platelet is represented as an IB object, i.e., as a 
collection of elastically-linked Lagrangian points that each move at the local fluid velocity. 
New elastic links are created dynamically to model the adhesion of a platelet to the injured 
wall or the cohesion of activated platelets to one another. Multiple links can form between 
a pair of activated model platelets or between a model platelet and the injured wall, and 
these links collectively represent the ensemble of molecular bridges binding real platelets to 
one another or to the damaged vessel. The links exert forces on the surrounding fluid to 
resist motions which would otherwise separate the linked entities. Links may break if subject 
to sufficiently high stress. Model variables are fully coupled: the fluid carries the activator 
and platelets, while the interplatelet forces, potentiated by chemically-induced activation of 
the platelets, determine the local flow. In this paper, we focus on mechanical interactions, 
not the activation process, and so we specify the conditions under which a platelet becomes 
activated and able to cohere with other activated platelets. 

We model platelets as closed curves of interconnected IB points in 2D and closed surfaces 
of interconnected IB points in 3D. A platelet's area or volume is determined by the region 
enclosed by the curve or surface and is preserved because of the incompressibility of the 
fluid. We describe these details for both the traditional IB method and our proposed new 
IB method in later sections. 



1.3 Motivation for the RBF-IB method 



We are motivated by the application of the IB method to platelet aggregation in blood 
clotting. Real platelets circulate with the blood in an inactive state in which they have a 
discoidal shape. In IB modeling, inactive platelets are approximately elliptical or ellipsoidal 
in 2D and 3D, respectively, while activated platelets are approximately circular in 2D and 
spherical in 3D. Piecewise linear approximations of platelets are currently used in IB methods 



applied to the simulation of platelet aggregation {e.g. [5,9,10]). 

In [24], we explored alternative methods for the modeling of platelets, focusing on methods 
that decreased the computational time necessary to maintain and update platelet geome- 
try and motion and had comparable or better error characteristics to the standard models. 
Specifically, we compared interpolation with Fourier-based techniques (trigonometric poly- 
nomials in 2D and spherical harmonics in 3D) and interpolation with radial basis functions 
(restricted to the unit circle in 2D and unit sphere in 3D). We found that interpolation with 
radial basis functions offered accuracy and computational cost comparable to that offered 
by Fourier-based methods (and better than that offered by standard finite-differences and 
piecewise-quadratic interpolation used in conjunction with traditional platelet simulations) 
in modeling a target shape, its normals and tension forces computed on its surface when the 
target shape was infinitely smooth. Furthermore, interpolation with radial basis functions re- 
sulted in better accuracy than that offered by both Fourier-based methods and the standard 
combination of finite-differences and piecewise-quadratic interpolation used in many versions 
of the IB method when the target shape had only one or two underlying derivatives. In this 
situation, use of radial basis functions led to a computational cost comparable to that of 
Fourier-based methods and lower than that of the standard combination of techniques used 
in many IB methods. 

In [24], we established that interpolation with radial basis functions is superior to popular 
alternatives in analytic test cases not involving any fluid-platelet interaction. We now turn 
our attention to exploring the consequences of using a parametric RBF geometric model 
within the full IB method, with platelet aggregation as our target application. In this work, 
we propose a new immersed boundary algorithm that utilizes the features afforded by our 
RBF geometric model. 

For results, we will compare the RBF-IB method to the traditional IB method for shapes 
that typify observed platelet geometries. We quantify the differences between the RBF-IB 
and traditional IB methods by comparing positions of platelets as generated from these two 
methods. We also determine the extent of volume loss and the maximum time-step size 
allowed by each method. We then provide timing comparisons between the two methods as 
a function of grid size, the number of immersed boundary points on each platelet's boundary 
and the number of platelets in a simulation. Our results show that the RBF-IB algorithm 
offers significant savings in terms of computational cost without a significant difference (in 
terms of platelet positions) from the traditional IB method. We furthermore find that the 
RBF-IB method allows for larger time-step sizes and lower volume loss than the traditional 
IB method, even with the very high stiffnesses required in platelet simulations. We then 
present results from a 2D platelet aggregation simulation computed using the RBF-IB and 
the standard IB methods. Finally, we present a simple example in 3D of two platelets as 
simulated by the RBF-IB method interacting with the flow, with each other, and with a 
wall. 

The paper is organized as follows. In Section 2 we briefly discuss the traditional Immersed 
Boundary Method for simulating fluid-structure interaction. In Section 3 we review the dif- 



ferent geometric modeling strategies that we seek to compare: piecewise linear and RBFs. In 
Section 4 we review the components necessary for handling immersed elastic structures in 
the IB method. In Section 5 we provide implementation details, both in terms of computing 
platelet elasticity and in using the constitutive models within an Immersed Boundary algo- 
rithm. In Section 6 we present our comparison of the RBF-IB method with the traditional 
IB method. We also provide energy estimates for RBF-IB simulations. We then present re- 
sults from a large platelet aggregation simulation in 2D and a simple platelet aggregation 
simulation in 3D. Section 7 contains a summary of our findings and a discussion of future 
research directions. 

Notation : We denote vectors with as many components as the spatial dimension (two or 
three) in bold. We denote vectors with as many components as the number of data sites 
(A'd) or sample sites (A's) by underlining. We indicate matrices with (A^d) or (A's) rows and 
two or three columns in bold with underlining. 



2 Review of the Immersed Boundary Method 



The IB method was developed to study the interactions of a viscous incompressible fluid 
with one or more moving and/or deformable elastic objects in contact with that fluid. The 
motion of the fluid influences the motion of the elastic objects and vice versa, and so the IB 
method involves coupled equations of motion for both types of material (fluid and elastic) and 
solves for both motions simultaneously. To review the IB method, we focus on a simple two- 
dimensional model problem in which a single fluid-filled closed elastic membrane is immersed 
in a viscous fluid. 

The physics of the model problem is that an elastic membrane is under tension and exerts 
forces on the adjacent fluid. These forces may cause the fluid to move and, correspondingly, 
cause points on the membrane to move along with the fluid. In the IB method, the fluid 
is described in the Eulerian frame through a velocity field u{x, t) and pressure field p{x, t) 
defined at every point x in the physical domain Q. The elastic membrane is described in the 
Lagrangian frame. Let the elastic membrane be parameterized by g, and denote by X{q,t) 
the spatial coordinates at time t of the membrane point labeled by q. The IB equations are 
the following coupled equations of motion for the fluid variables u{x,t) and p{x,t) and the 
membrane configuration X{q,t). 

p{ut + u-Vu) = -Vp + fiAu + f, V-w = 0, (1) 

F{q,t) = F(x{q,t),-^X{q,t)Y (2) 

fix, t) = j F{q, t) 6{x - X{q, t)) dq, (3) 

dX r 

—-{q, t) = / u{x, t) 6{x - X{q, t))dx. (4) 

ot Jn 



Equations (1) are the Navier Stokes equations which describe the dynamics of a viscous 
incompressible fluid, of constant density p and constant viscosity /i, driven by a force density 
/ which here arises because of the elastic deformation of the immersed membrane. Equation 
(2) specifies the elastic force (per unit q) at each point of the immersed boundary object. The 
functional dependence of this force on the state of the boundary is specified appropriately to 
the material being modeled. An example is given below. Equation (3) defines the fluid force 
density f{x,t) in terms of the immersed boundary elastic force density F. By integrating 
both sides of this equation over an arbitrary region of the fluid, we see that the total fluid 
force on this region equals the total elastic force along the portions of the immersed boundary, 
if any, that pass through this region. The fluid force density is therefore concentrated along 
the immersed boundary curve. Equation (4) specifies that the velocity of each immersed 
boundary point equals the fluid velocity at the same location. This is a formulation of the 
no-slip boundary condition for viscous flows. The key idea in this formulation that makes the 
IB approach so useful in modeling biofluid problems is that as far as the fluid is concerned, the 
immersed objects are seen only through the force field /. Even if the objects move or deform 
substantially, there is no change in the geometry of the fluid region; fluid is everywhere and 
only the distribution of forces exerted on the fluid by the elastic objects changes. In the 
model problem and the platelet applications, we assume that the IB objects are neutrally 
buoyant; the IB membrane itself carries no mass each object's mass is attributed to the fluid 
in which it sits. 

For this study, we restrict our attention to the Forward Euler update scheme for both the 
traditional IB method and our proposed new method, since most traditional IB methods use 
an explicit integration scheme to update IB point positions. 



3 Geometric modeling of platelets 



In this section, we review the two geometric modeling strategies to be compared in the context 
of the Immersed Boundary method applied to platelet aggregation. For a full description of 
these strategies, see [24]. 



3. 1 Piecewise linear model 



In the traditional IB method, the surfaces of platelets are represented by a collection of 
Immersed Boundary points. We henceforth refer to the IB points in the traditional IB method 
as sample sites, and denote them by -X^s(^) = ^{(1, t) for g G F and at a particular time t. The 
surface elastic forces of the platelets are spread from these sample sites into the neighboring 
fluid. Both tension and bending forces are computed using a collection of springs (typically 
linear springs) between pairs of sample sites. An explicit piecewise linear interpolant of 
the surface is not formed. If other information (such as normal vectors) is needed at the 
sample points, an approximation to the surface represented by the sample sites may be 



formed by a piecewise quadratic interpolation of the sample sites {e.g., [29]). After the 
incompressible Navier Stokes equations are solved, velocities from the portions of the Eulerian 
grid surrounding the sample sites are interpolated to the sample sites and used to move the 
platelets. 



3.2 Parametric RBF model 



The RBF method is a popular tool for approximating multidimensional scattered data. An 
excellent overview of the theory and application of this method can be found in the two books 
by Fasshauer [3] and Wendland [26]. The restriction of the RBF method to interpolation on 
a circle and or sphere began to receive considerable attention from a theoretical standpoint 
in the mid 1990s (see [4, §6] for a discussion). When restricted to these domains, the RBF 
method is sometimes referred to as the zonal basis function (ZBF) or spherical basis function 
(SBF) method in the literature [26, Ch. 17]. We however use the more popular term RBF to 
describe the interpolation technique. Several studies have provided error estimates for RBF 
interpolation on circles and spheres; see, for example, [18,19]. In [18], it was shown that 
these interpolants can provide spectral accuracy provided the underlying target function 
is sufficiently smooth, while in [19], error estimates are given in the case that the target 
function belongs to a Sobolev space. The RBF method has also been used successfully for 
numerically solving partial differential equations on the surface of the sphere [6,7], as well 
as more general surfaces [11,22]. 

Here, we present the RBF model developed in our earlier work [24]. It is based on explicit 
parametric representations of the objects. Since our target objects are platelets, which in 
2D models are nearly elliptical or circular (ellipsoidal or spherical in 3D), we choose a polar 
parameterization. We use our model to define operators necessary for the computation of 
geometric and mechanical quantities required by the IB method. 



3.2.1 Two-dimensional platelet modeling 

In 2D, we represent a platelet surface at any time t parametrically by 

X(A,t) = (X(A,t),y(A,t)) (5) 

where < A < 27r is the parametric variable and X(0,t) = X(27r,t). We explicitly track 
a finite set of N^ points X^{t), . . . , Jfgl^)) which we refer to as data sites. Here X At) := 
X(A^, t), j = 1, . . . , A'd, and we refer to the parametric coordinates A'j', . . . , A^^ as the data 
site nodes (or simply nodes). We construct each component of X by using a smooth RBF 
interpolant of the data sites in parameter space as discussed in detail below. 

We also make use of derivatives of the interpolant at the data sites and we use the inter- 
polant and its derivatives at another set of prescribed sample points or sample sites, which 
correspond to N^ parameter values: A^, ..., A^ . Here we show how to form the interpolant. 



how to evaluate its derivatives at data sites, and how to evaluate the interpolant and its 
derivatives at sample sites. Note that the set of sample site parameter values may be disjoint 
from the set of data site ones. 

Here, we explain how to construct an RBF interpolant to the X component of X using the 
data {Xf,Xf(t)), ..., (A^ ,X^ (t)); the construction of the Y component follows in a similar 
manner. Let 0(r) be a scalar-valued radial kernel, whose choice we discuss below. Define 
X(A,t)by 



X(A,t) 






k=l 



2-2cos(A- A^; 



(6) 



Note that the square root term in Equation (6) is the Euclidean distance between the points 
on the unit circle whose angular coordinates are A and A^. For this paper, we consider two 
families of radial kernel functions, the multiquadric (MQ) for 2D platelets and the inverse 
multiquadric (IMQ) for 3D platelets. These are given explicitly by 



MQ: 
IMQ: 



[r) 
(r) 



1 + (er) 



\/i + i^^y 



(7) 
(8) 



where e is called the shape parameter. To have X{X, t) interpolate the given data, we require 



that the coefficients c^, k 



Nd satisfy the following system of equations: 



(^JVd.i 






(rN,, 



Nd) 





'x' 




X,d(t)" 




^2 


= 


Xf(t) 








Xl(t)_ 



Cd 



(9) 



lUt) 



where r^-^fc = w2 — 2 cos(A^ — Af). Since r^^k = fkj, the matrix A in this system is symmetric. 
More importantly, for the MQ and IMQ kernels, A is non-singular (and positive-definite for 
the IMQ) [3,26]. 

In our application, we want to be able to evaluate X(A,t) at sample sites corresponding to 
parameter values A^, ..., A^_^, that stay fixed over time. While we could use Equation (6) to do 
this, it is much more convenient from a notational and computational perspective to construct 
an evaluation matrix that combines the linear operations of constructing the interpolant to 



X^(t) = [2Ld{t)i Y^{t)], for any t, and evaluating it at A^, 






The evaluation matrix can 



be constructed by first noting that Equation (6) can be written as 



X{X,t) 



('y^2-2cos(A-Af)') ■ ■ ■ (y^2 - 2cos(A - A^j) 



cl 



(10) 



6(A) 



T 



Since c^ = A ^2Ld{'t)j we can further write Equation (6) independent of c^ as X(A,t) 
b{X)^ A~^ 2Ld{'t) ■ The evaluation of X(A,t) at A^, ..., A^ can then be obtained as follows: 



X(A^t) 




'biXlf 








= 




A-^X^{t) - 


= BA_^2U{t) 


X{X%^,t)_ 




b{xkr_ 




Ss 


XJt) 




B 







;ii) 



So, given the data sites 2Ld{'t) ^t any time t, we can interpolate these values with an RBF 
expansion and evaluate the interpolant at the sample site nodes A^, ..., A^_, to get Xg(t) by 
the matrix-vector product £^X^(t). In fact, this same procedure can be used to give values 
at sample site nodes for any quantity whose values we have at data site nodes and which 
we represent using an RBF expansion {e.g., Y_^(t) = [Yf(t)---Y^^{t)]^). Furthermore, the 
evaluation matrix S^ can be precomputed once at t = and used for all subsequent times. 



We also need to compute geometric quantities such as tangent vectors, and mechanical 
quantities such as forces at data sites and/or sample sites. These quantities require computing 
derivatives with respect to A of the platelet surface {X{\, t), Y{X, t)). We use the RBF-based 
representation of the surface to compute these derivatives. As in the case of interpolation 
and evaluation, we will express derivatives of the RBF interpolant in matrix-vector form. 
Toward this end, we use similar notation to Equation (10) and define the vector 



^a(A) := ^hiX) 



QV 



A=A 



(9A" 



2-2cos(A-Af) 



A=A 






2-2cos(A-Ai ) 



A=A. 



for some < A < 27r. Just as b{X)^A ^2Ld{'t) gives the value of X{X,t), we can use 6"(A) to 
obtain the n^^ derivative of X{X,t) with respect to A as 



Qn 

dXj^ 



X(A,t) 



mxyA-'2Ld{t). 



A=A 



The evaluation of the n^^ derivative of X{X, t) at the data site nodes A^, . . . , A^ can then 



be obtained as follows: 

d 



Qn 



X{X,t) 



9A" 



X(A,t) 



A=A"j' 



A — A»r 



n/ \d\T 



hliXi 



m\%,r 



A-^2Ld{t) = B^,A-^Xd{t). 



(12) 



^Ad 



^Ad 



In a similar manner, the evaluation of the n^^ derivative of X{X, t) at the sample site nodes 
A^, . . . , A^ can be obtained by 



gn 

dX^ 



X{X,t) 



QU 

dX' 



-X{X,t) 



A=A5 



A=A= 





A-'X^) - 




Bl 







(13) 



For given data sites X^{t) at any time t, we can interpolate these values with an RBF 
expansion and evaluate the n*^ derivative of the interpolant at the data site nodes by the 
matrix- vector product V^^X^it) and at the sample site nodes by V^^Xdit). We refer to the 
^d X ^d matrices "D^d and the N^ x A'^d matrices "D^s as RBF differentiation matrices. 

As in the case of the evaluation matrix E^ defined in Equation (11), the matrices "D^d and V^s 
can be used to give values at respective data site or sample site nodes of the n^^ derivative 
of the RBF interpolant of any quantity whose values we have at the data site nodes {e.g., 
Y^{t) = [Y{^{t) ■ ■ ■Y^^it)]'^). Furthermore, these matrices can also be precomputed once at 
t = and used for all subsequent times. 



3.2.2 Three-dimensional platelet modeling 

In 3D, we represent a platelet surface using the following spherical parametric notation: 



XiX,e,t) = (X(A,^,t),F(A,^,t),Z(A,^,t)), 



(14) 



where —n < A < vr and — 7r/2 < 9 < 7r/2. Here, the end conditions on X in A are 
X(— TT, 6,ti) = X(7r, 6*, tj), while the end conditions in 6 are X(A, 7r/2, tj) = X(A + 7r, 7r/2, tj) 
for -TT < A < and X(A, -7r/2, t^) = X{X - vr, -7r/2,ti) for < A < tt. Similar to 2D, 
we track a finite set of A^d data sites Xi{t), . . . ,X2{t), where we now define Xj{t) : = 
X(A^, 9j,t), j = 1, . . . , Ad and refer to the parameter pairs (A^, 9j) as data site nodes. We 
again construct each component of X by using a smooth RBF interpolant of the data sites in 
parameter space. For the spherical parameter space the RBF interpolant of the X component 



10 



of X is given by 



iVd 



x{\e,t) = Y^ 



ji . 



k=l 



2(1 - COS e COS ^^ cos(A - A^) - sin 6 sin O"^] 



(15) 



where is defined as in the 2D platelet model. The argument of (p is the Euchdean distance 
between the points on the unit sphere with spherical coordinates (A,^) and (A^,^^). The 
coefficients {c^} in Equation (15) are found by solving a similar linear system to Equation 
(9): 



AcJ 



2U{t), 



(16) 



where c^ is a vector containing the interpolation coefficients and 2Ld{'t) is vector containing 
the X data sites. The entries of the matrix A in Equation (16) are given by Aj^^ = 4'{rj,k), 
with Tj fc := 



'2(1 — cos 6*^ cos 6^^ cos (A^ — A^) — sin 6*^ sin 6^^). This matrix is also guaranteed 
to be non-singular for the MQ and IMQ kernels (and positive definite in the case of IMQ). 
The Y and Z components of X are interpolated in a similar manner. 

As in the 2D model, we want to construct an evaluation matrix that can be used to form an 
RBF interpolant of any quantity whose values we have at the data site nodes and evaluate it 
at a collection of A^'s sample site nodes (A^, 91), ... , (A^^, 0%J. Following a similar procedure 
to the 2D model, we define a vector function bx{X, 0) with entries 



ib{X,e))j := ( y^2(l - cos e cos Oj cos(A - A^^) - sin 6 sin Of) j , j 



1 N.. 



(17) 



Given the data sites Xjj(t), we can use b{\,9) to write Equation (15) as X{\,9,t) = 
b(X,6)^A~^2Ldi't)- The evaluation of the interpolant of Aj(t) at the sample site nodes is 
then given by 



X{Xl,9l,t) 



s\T 



kixi^ei 



h.{X% ,9lfV 



A~^2Ldit) = BA-^2Ld{t)- 



2L{t) B 

The iVg X N(i matrix Sg is the evaluation matrix for the 3D platelet model. 

RBF differentiation matrices can also be constructed for the 3D platelet model in a similar 
fashion to the 2D model for computing (mixed) partial derivatives of the RBF interpolant 

Qm+n 

of some quantity defined at the data site nodes, e.g., ^^^ ^^^ X{X, 9, t). Toward this end, we 



dX"'d9' 



define a vector function 6^g"(A, 9) with entries 

am+n 



{bTfiKO))j 



QXmQQn 



ibiX,9)), 



J=l, 



m 



X=\j 



11 



where b{X,e) is defined in Equation (17), -tt < A < vr and -7r/2 < 6* < 7r/2. The A^d x N^ 
RBF differentiation matrices for evaluating mixed partial derivatives at the data site nodes 
are given by 

^\dgd ■— J->xdQd^ , {^■Jj 

where the i^^ row of B^/'g^ is equal to 6^g"(Af , OfY, i = 1, . . . , A^d- Similarly, the A^s x ^d 
matrices for evaluating mixed partial derivatives at the sample site nodes are given by 

where the i^^ row of -B^'^s is equal to llx'ei,^,, ^1)"^; ^ = 1, • • • , ^s- 

As in the 2D platelet model, the evaluation matrix E^, and differentiation matrices ^^'gd 
and P^'gs only need to be computed once at t = and stored, and can then be used for all 
subsequent times. 



4 Lagrangian force modeling 



In this section, we present the components needed for modeling Lagrangian forces, i.e., 
the forces on the surface of platelets due to their individual configurations, along with the 
tension and bending models that we utilize in our implementation of the traditional IB and 
the RBF-IB methods. 



4.1 Components for the force model 



We represent a 2D platelet parametrically by X{\,t) as discussed in the previous section 
and define 

T := g^XiX,t) = (—X{\,t),—Y{\,t)\ = {Tx,Ty). (21) 



The unit tangent vector to X{X,t) is then given by 

'r- = 7^ = irx,ry). (22) 

||t|| 

For a 3D platelet, we use the parametric representation X{X,9,t) discussed in the previous 
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section and define 

t' --—xix, e, t) = (-^x{x, e, t), |-y (a, e, t), ^z{x, e, t)] , (23) 

The unit tangent vectors to X{X, 6, t) are then given by 

T^:=-— and f'^ := ^, (25) 



while the unit normal vector is given by 



^ ■= T^J^V (26) 



4-2 Force model in 2D 



In our experiments, we assume that the Lagrangian force at a point on a platelet is the sum 
of a tension force, a bending-resistant force and possibly a force due to a bond between that 
point and a point on another platelet or the vessel wall. For the tension force, we use the 
fiber model defined in [21], according to which the elastic force density at X(Ai,tfc) due to 
tension is given by 



F^{\M)=§^iTT) 



(27) 



where T = fct(||T|| — ^o) is the fiber tension, where /ct > is constant. We set /q.j = II '^11 1 a, t ' 
where to is the initial time of the simulation. 

For a bending-resistant force, we use a linear variant of the force defined in [13] and define 
the elastic force density at X{Xi,tk) due to how much the platelet surface there is bent to 
be 



F^{KM) = -kJ^'^ ^'^°^ 



(28) 



^ dX^ dX^ 

Here X" = X(Aj,to) is the initial configuration of the platelet and /cb > is constant. 

We defer discussion of how we compute the forces given by Equations (27) and (28) to 
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the next section, since the implementation is different for the RBF and piecewise-hnear 
representations of the platelet boundary. 

However, the force acting on a platelet due to other platelets (and/or walls) is common 
to both methods. We use the spring force defined in [10]: let Pi,P2, ---yPNp be the indices 
corresponding to the platelets in the domain. Let pi and p2 be the indices of two platelets 
which are linked at sample sites Xp^(AU and Xp^{X^J. The force at Xp^(AU is then given 
by: 

where Kq and /o,c are the interplatelet cohesion spring stiffness and the resting length, 
respectively. We use asimilar formulation for platelet-wall links. 



4.3 Force Model m 3D 



As in the 2D force model, the total Lagrangian force in 3D is expressed as the sum of a tension 
force, a bending-resistant force, and possibly a link force. Since platelet surface membranes 
(like all cell membranes) are lipid bilayers (see [1]), we use tension and bending models, based 
on curvaturem developed for lipid bilayer membranes as in [2], with modifications introduced 
to take into account the existence of an initial shape. 

The force due to surface tension is given by 

F^{\, e,, tk) = -7 {'2Hf]\,^ - 2HfjQ , (30) 

where 7 > is the coefficient of surface tension. 

The bending-resistant force is given by 

F^(A„ e,, tk) = h (2H{H^ -K)f]^ - 2H{H^ - K)f] ] , (31) 






to 



where /cb > is the bending rigidity. H is the mean curvature of the surface and K is the 
Gaussian curvature [12, §16.5]. H and K may be calculated as follows: 

_ eG - 2fF + gE 
^~ 2{EG-F^) ' ^^^' 



and 



•'-wr&^' (33) 
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where E, F, and G are coefficients of the ffist fundamental form, 

E = T^.r^^ F = t^-t', G = t' ■ T^ (34) 

and e, /, and g are coefficients of the second fundamental form, 

' - (H ■^'•f- (H ■ "• « - (m^') ■ ^ '^=> 

The inter- platelet force in 3D is given by the 3D analog of Equation (29). 
5 Implementation 



In this section we present the implementation details for evaluating positions on the La- 
grangian objects and computing their Lagrangian forces as presented in the previous sec- 
tion. For the piecewise linear representation, these forces are computed at the IB points. For 
the parametric representations using the RBF model, these quantities are computed at the 
sample sites, which do not necessarily correspond to the data sites. Finally, we present the 
RBF-IB algorithm that incorporates our RBF model for platelet elasticity. 



5. 1 Piecewise linear models 



Traditionally, finite-difference approximations of Equations (27) and (28) are used in con- 
junction with piecewise linear methods in 2D {e.g. [13]). We use a second-order central 
difference to discretize the derivatives involved in the computation of both the tension and 
bending forces (including tangent lengths). It is useful to think of these finite difference 
approximations to the constitutive model as Hookean springs connecting pairs of IB points. 



5.1.1 The 2D piecewise linear Immersed Boundary algorithm 

First, we present the steps of the traditional Immersed Boundary algorithm. For a more 
detailed discussion, see [21]. 

(1) The resultant Fq of all of the force contributions that act on an IB point Xq is calculated 
for each q. 

(2) These forces are distributed to the Eulerian grid used for the fluid dynamics equations 
using a discrete version of Equation (3): 

f{x,) = J2 FM^, - Xq)dq. (36) 
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Here, Xg and Xg are the coordinates of grid point g and IB point q, respectively, Fg 
is the Lagrangian force (per unit q) on the IB point, dq is the increment in parameter 
q between consecutive discrete sample sites, and 6h is a discrete approximation to a 
two-dimensional (5-function. 

(3) With the fluid force density fg now known at each grid point, the fluid velocity is 
updated taking one step with a discrete Navier-Stokes solver. We use a projection- 
method [15]. 

(4) Denoting the new velocity field by u'^'^^, IB points are updated by a discrete analog of 
Equation (4) 

X^°- = Xg + AtUg = X, + AtY^ ur^hixg - X,)h\ (37) 

9 

where h is the fluid grid spacing and 5h is the same approximate (5-function as used in 
Equation (36). 

For the approximate (5-function, we use the "cosine" form described by Peskin [21] which 
ensures that the entire IB force is transmitted to the grid, that the force density on the grid 
is a continuous function of the IB point locations, and that the communication between grid 
and IB points is very localized. 

After the above steps are completed, new links are formed and existing ones are broken using 
the model's rules for these types of events. 



5.2 RBF models 



We now turn to a discussion of the changes introduced to the traditional algorithm by 
incorporating an RBF model. We first explain how we select node sets to construct the RBF 
operators presented earlier; then, we explain how to implement these operators in an efficient 
fashion. Finally, we present our RBF-based algorithm for computing platelet forces. 



5.2.1 Node and evaluation points 

In order to construct the operators utilized by our algorithm, we must first choose an appro- 
priate node set. In 2D, we use A^d equally-spaced points on the interval (0, 2n] as the data site 
node set {)^'k]k=v This gives a uniform sampling in the parametric space. More importantly, 
since our target objects are either circular or elliptical in 2D, these node points correspond 
to a nearly uniform distribution of data sites on the object. We also use As >> Ad equally- 
spaced points in the interval (0, 27r] as the set of sample site nodes {Ajl^Ii since this results 
in a set of sample sites that are well distributed over the object. 

To get a good sampling of our nearly ellipsoidal or spherical objects in 3D, we don't use 
equally spaced points in the spherical coordinate system as our data site nodes {(A^, 6^^ )}^^^ 
because they result in a non-uniform clustering of data sites at the "poles" of the spherical 
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coordinate system. Instead we use node sets that give a quasi-uniform distribution of data 
sites on the unit sphere. Since only a maximum of 20 points can be evenly distributed on a 
sphere, there are a myriad of methods to define and generate a quasi-uniform distribution 
for larger numbers of points [16]. We use so called minimal energy (ME) point distribution 
for our RBF models discussed in [27] and freely available for downloading of various sizes 
from [28]. The ME points are generated by finding a distribution of nodes that minimize 
an electrostatic type energy potential. For the set of sample site nodes, {(A^-, O^j)}^!^, we use 
A's >> Nd_ ME points, which again results in a well distributed set of sample sites on the 
object. 



5.2.2 Evaluation and differentiation operators 

We have formulated our operators to ensure that operations like evaluation of the interpolant 
and computing derivatives (and therefore the constitutive model) do not require solving 
a linear system for any time step of the platelet simulation except the initial step. This 
is possible because, though the data sites and sample sites move over the course of the 
simulation, their values in parameter space do not change. 

For the 2D RBF model of the platelets, the evaluation matrix S^ in Equation (11) and dif- 
ferentiation matrices V^^ and V^s in Equations (12) and (13), respectively can be computed 
using the FFT as discussed in [24]. This is possible since the data site nodes {A^l^^^ are 
equally spaced, which results in the A matrix defined Equation (9) having a circulant ma- 
trix structure. For the 3D RBF model of the platelets the interpolation matrix A defined 
in Equation (16) is not circulant. However, since we use the IMQ kernel, A is symmetric 
positive definite. Thus, we compute the Cholesky factorization of A once, and use it to com- 
pute the the evaluation matrix Ss in Equation (18) and differentiation matrices ^^'^d and 
D^'gs in Equations (19) and (20), respectively. The costs and accuracy of the RBF models 
are elaborated upon in the discussion of the 2D results. 



5.2.3 Algorithms for computing platelet forces with RBFs 

We now describe the implementation of the constitutive models outlined in Section 4. We 
present algorithms for computing platelet forces in both 2D and 3D. 

Notation : In the description of the algorithms below we use standard matrix-vector op- 
erations such as multiplication as well as non-standard operations like element-by-element 
multiplication of matrices and vectors (sometimes called the Hadamard product). We denote 
this latter operation with the o operator. For example, if J and L are A^d x 2 matrices and 
^ is a vector of length A'd then the i^^ row oi J_o L_ and RoJ_ are given by 

(JoX),,i;2= [(J)^,l(^)^,l, (J).,2(^).,2] 

where i = 1, . . . , N^. 
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We define t_^ = Vj^^Xj^it) , the A'j x 2 matrix of tangent vectors at the data sites at time 
t and ||Td||, the N^ vector containing the two-norm of each row of t^. The algorithm for 
computing platelet elasticity in 2D is as follows: 

(1) Initialization (t = to): After creating and storing the RBF evaluation matrix as in 
Equation (11) and differentiation matrices as in Equations (12) and (13), compute for 
each platelet: 

(a) The rest lengths for the tension force at the data sites: /o = Zd = 'Dy, X, j(tn). 

(b) The bending-resistant force term for the platelet 's initial configuration at the data 
sites: Vj^sX^Un). 

(2) For each time step {t = tj^, k > 1), compute for each platelet: 

(a) The length of the tangent vectors r^ = Vl^X^Uk) at the data sites: ||Td||; and the 
unit tangents at the data sites: r^- 

(b) The tension at the data sites, using the constitutive model: T^ = /^tdlTdll — ^)- 

(c) The tension force at sample sites: £g = V\sZ_^, where Z_^ = Td° i-d- 

(d) The bending force at sample sites: F^ = —k\, {V\sX^{th) — 'D\sX^{tn)). 

(e) The interplatelet cohesion force from Equation (29) at the sample sites: F_^ . 

(f ) The total Lagrangian force at the sample sites: F_^ = FJ + F^ + F_^. 

To compute platelet forces in 3D, we require an algorithm to compute the mean curvatures 
ifg and the Gaussian curvatures K^ at the sample sites at any given time. The algorithm is 
as follows: 

(1) Compute the tangent vectors at the sample sites, r^ = V^s ^sX^it) and rf = V^^, g^X^it). 

(2) Compute the unit normal vectors at the sample sites, r) , by computing the cross prod- 
ucts of each row of r^ with the corresponding row of r^ and then normalizing the length 
of the resulting vector. 

(3) Compute the coefficients of the first fundamental form at the sample sites, ^g, £s! ^^^ 
Gg according to the equations in (34). 

(4) Compute the coefficients of the second fundamental form at the sample sites, Cg = 
(I^A^.eXd(t)) o^, /^ = (p°L%Xd(t)) o^^, and^^ = (PlLV^d(t)) °i,- 

(5) Compute H_s and K_^, the curvatures at the sample sites, using Equations (32) and (33), 
respectively. 

With this algorithm in hand, we now present the algorithm for computing platelet elasticity 
in 3D: 

(1) Initialization (t = to): After creating and storing the RBF evaluation matrix as in 
Equation (18) and differentiation matrices as in Equations (19) and (20), compute for 
each platelet: 
(a) The initial mean curvature at the sample sites: H^. 
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(b) The initial mean Gaussian curvature at the sample sites: K^ . 

(c) The surface tension force terms at the sample sites for the platelet 's initial config- 
uration: — 7(2i7g) o r) . 

(d) The bending-resistant force terms at the sample sites for the platelet's initial con- 
figuration: /cb (2i7g o [Hj, o Hj, — Kj,)) o r) . 

(2) For each time step {t = tk, k > 1), compute for each platelet: 

(a) The mean curvature at the sample sites: H^ . 

(b) The Gaussian curvature at the sample sites: K^ . 

(c) The tension force from Equation (30) at the sample sites: F_^ . 

(d) The bending force from Equation (31) at the sample sites: F^ 

(e) The interplatelet force from the 3D analog of Equation (29) at the sample sites: 

(f) The total Lagrangian force at the sample sites: F^ = FJ + F^ + F_^. 



5.2.4 The RBF Immersed Boundary algorithm 

Equipped with the ability to evaluate constitutive models with our parametric RBF model 
of platelet surfaces, we now present the final Immersed Boundary algorithm, incorporating 
the changes introduced by the RBF paradigm. We note that we use the previously presented 
notation when operating on A'^^ x 2 or A^s x 3 matrices of sample sites {e.g., X^)- We refer to 
individual sample sites with index q as {Xs)q and individual data sites with index j as (-X'rf)^ 
when appropriate. This helps when comparing the RBF-IB algorithm with the traditional 
IB algorithm. 

(1) The total force on a platelet F^ is calculated using the appropriate algorithm from the 
previous subsection. 

(2) These forces are distributed to the Eulerian grid used for the fluid dynamics equations 
using a discrete version of Equation (3): 

f{x,) = J2FM^s-iX,),)dq. (38) 

Here, Xg and {Xs)q are the coordinates of grid point g and sample site q, respectively, 
Fq is the Lagrangian force (per unit q) on this point, dq is the increment in parameter 
q between consecutive discrete sample sites, and 6h is a discrete approximation to a 
(5-function. 

(3) With the fluid force density fg now known at each grid point, the fluid velocity is up- 
dated taking one step with a discrete Navier-Stokes solver. We use the same projection- 
method type of Navier-Stokes solver as before. 

(4) Denoting the new velocity field by it°°™, update the data sites {X^jj by a discrete 
analog of Equation (4) 

(^r ). = (^d), + At(t/d), ^ (Xd), + At ^ u^^-6,ixg - {X^),)h\ (39) 

a 

(5) Generate a new set of sample sites ^s('^ncw) by applying the RBF evaluation operator 
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to the data sites XT^ := Xi(^r, 



I.e., 



2L{tn 



^s2Ld [tn 



(40) 



In the results section, we explore the differences that arise from the traditional IB method 
when using this proposed algorithm. 



6 Results 



In this section, we present the results of our numerical experiments comparing the RBF-IB 
method and the traditional IB method. We express differences between the two methods by 
computing differences in sample point positions between objects (having started both simu- 
lations with the same set of sample points), centers of mass and timings. For brevity, we refer 
to the traditional IB method as the PL-IB method (Piecewise Linear-Immersed Boundary 
method). We note that since the PL-IB method does not represent a "gold standard" but 
simply the current working model, it is difficult to establish a notion of accuracy or calculate 
errors for the RBF-IB method. We first examine the behavior of both methods on simple 
test cases, select the shape parameter based on these tests and then run simulations with 
many platelets. We present these experiments and results in 2D, and after thoroughly study- 
ing our algorithm's behavior in 2D, show the results of a platelet aggregation simulation in 
2D. We then show the results of a 3D simulation involving platelets. Our rationale for not 
performing comparisons in 3D is that with the traditional IB method in 3D, constitutive 
models are rarely used; on the other hand, within the RBF-IB method in 3D, we have cho- 
sen to model the platelet surface as a lipid bilayer membrane. Comparisons between the two 
methods of the sort made for the 2D therefore cannot be made in the 3D case, since both 
the computational cost and the dynamics are different. 

All tests employ periodic boundary conditions on the left and right ends of the domain and 
no slip boundary conditions on the top and bottom. 

Finally, in the results below we let X^ and X^ denote the A^s sample sites of the RBF- 
IB and PL-IB methods, 
following two measures: 



respectively. Differences in these values are computed using the 



Normalized RMS difference: 
Normalized maximum difference: 



\ 



A^B 



A^. 



U(x 



RBF\ 



ix 



s Jj 



i=i 



max 



X 



KBF\ 



X% 



/Cp, 
/Cp, 



(41) 
(42) 



where Cp is the circumference of the object. Both are normalized measures of the distances 
between corresponding sample sites from the RBF-IP and PL-IB methods. In all the exper- 
iments, the sample sites for both the RBF-IB and PL-IB methods are initialized to be the 
same value. 
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Fig. 1. Difference in tlie reconstruction of tlie shape of the objects between the RBF-IB method 
and the PL-IB method for the Drift Test with a time step of At = 10~^ The figure shows the RMS 
normahzed difference in sample site positions according to Equation (41) at a time just before the 
objects leave the channel vs. the number of data sites A'd- The difference between the centers of 
mass of the objects was less than 0.1% of the channel height. The number of samples sites was set 
to Ns = 100. 

6.1 Verification tests 



In this section, we seek to verify the correctness of the RBF-IB method and quantify its 
differences from the PL-IB method. These tests were performed on solitary elastic objects 
that start as ellipses, with the same set of sample points for both methods. While neither of 
these cases are likely to be encountered in the context of platelet simulations, they provide 
simple and natural extensions of the tests performed on static platelets in [24] to incorporate 
interactions with the fluid. 



6.1.1 Drift test 

For this test, we initialize the simulations with the platelets on the inlet end of a channel. We 
then set up a forced flow with a parabolic velocity proflle. We disable the spreading of forces 
from the object to the fluid. As a consequence, the deformation of the objects is solely due to 
variations in the fluid velocity. The simulation was conducted until the elastic objects left the 
channel. At the time just before this exit, we measure normalized RMS differences in sample 
site positions of the solitary objects from each IB simulation according to Equation (41). 
Figure 1 shows these differences as the number of data sites N^ increases for two different 
grid sizes. We can see that the sample site positions from the RBF-IB method agree more and 
more with the positions from the PL-IB method with increasing A'^d. However, the differences 
between the sample site positions are larger on the finer grid than on the coarse grid. We 
note that this is not by itself a cause for concern since the PL-IB method does not represent 
a gold standard. 
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6.1.2 Time-step sizes and Area change on the Elasticity Test 




Fig. 2. The figure on the left shows the starting configuration of the elastic object, whose reference 
configuration is a circle. The figure on the right shows the same object after five units of simulation 
time. The number of datasites was set to A^d = 27 and the number of sample sites to Ng = 200. 
The figures were generated from the RBF-IB method. The figures generated by the PL-IB method 
are identical and hence are not shown. 

For these tests, we set the initial velocity of the fluid to zero. We set the reference configu- 
ration of an elastic object to be a circle, but the initial geometry to be an ellipse. We now 
allow the object to spread forces into the fluid. The object should therefore deform from 
an elliptical shape to a circular shape of equal area. The objects in these simulations are 
initially centered at (0.5,0.5). The test was conducted for 10^ time-steps on both a 32 x 64 
grid and a 64 x 128 grid. This is enough time-steps for the objects to reach steady state. 
The initial and final shapes of the object are shown in Figure 2. In this setting, we measure 
the percentage area change in the object as a function of time-step size for each method. 
We also use this experiment to measure the sensitivity of the RBF-IB method to the shape 
parameter e in the MQ kernel given in Equation (7). 
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Fig. 3. Percentage change in area in the RBF-IB method and the PL-IB method for the elasticity 
test for N([ = 27 and iVj = 54 as a function of the time-step (dt=At), measured after steady state 
has been achieved. The figure on the left shows the percent area loss on a 32 x 64 grid for both 
methods. The figure on the right shows the percent area loss on a 64 x 128 grid for both methods. 
The shape parameter was set to e = 3.1 for both N^ = 27 and N^ = 54. The number of sample 
sites was A^s = 200. 
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Figure 3 shows that with A'd = 27, the object simulated by the RBF-IB method has a greater 
change in area (after steady state has been achieved) compared to its initial configuration 
than when it is simulated by the PL- IB method. However, the area conservation of the RBF- 
IB method is better with N^ = 54 than the PL-IB method. This is a clear indication that 
there is a lower limit to the number of data sites one can use within the RBF-IB method, 
similar to the lower limit on the number of data sites established in [24] based on geometric 
accuracy. However, it is also clear from this test (the case of N^ = 54) that one can still use 
fewer data sites within the RBF-IB method than IB points in the traditional IB method and 
still see improvements over the latter. 

Although results are not presented, we have also found that these trends persist on finer 
grids, albeit with smaller percentage area change. We note that change in area is dependent 
on how well the grid resolves the object as well as the magnitude of the stiffnesses of the 
object. For both grid sizes, we used the same values of these parameters. When the stiffnesses 
are reduced down to the values presented in [14], the changes in area observed with N^ = 27 
and A^d = 54 both agree with the results presented there. 

These tests also show that the RBF-IB method allows stable simulations with time-steps 
that are significantly larger than those allowed by the PL-IB method. This is dependent on 
the number of data sites used in the RBF-IB method. We show the change in area as a 
function of the size of the time-step At for two values of A^^i in Figure 3. With the PL-IB 
method, we were unable to simulate with a time-step larger than At = 3 x 10~^ on the 
32 X 64 grid and At = 2 x 10"^ on the 64 x 128 grid. With the RBF-IB method, however, 
we were able to simulate with a time-step of even At = 6.5 x 10~^ on both these grids. 



6.2 Selecting the shape parameter for platelet simulations 
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Fig. 4. The normalized maximum difference in sample point positions of the objects from the 
RBF-IB and PL-IB methods vs the shape parameter, e. The differences are computed according to 
Equation (42). The figure on the left shows the results for the drift test and the figure on the right 
shows the results for the elasticity test. The number of datasites was set to N^ = 27. 
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In earlier work [24], we examined the impact of the RBF shape parameter e on shape recon- 
struction for a problem without fluid-structure interaction. We now study the effect of the 
shape parameter on differences between the RBF-IB method and the PL-IB method that 
could arise based on both motion and deformation. We measure the maximum differences 
in sample site positions according to Equation (42) as a function of the shape parameter in 
both the drift test and the elasticity test. The results from the tests are shown in Figure 4. 

It is easy to see from Figure 4 that the differences in sample site positions between the 
two methods are largely insensitive to the shape parameter, both in the drift test and the 
elasticity test. This is fortunate, as the selection of a shape parameter is generally a non- 
trivial task. The insensitivity of the RBF-IB method to the shape parameter is especially 
notable in the elasticity test, where the object undergoes more deformation than one expects 
to see in a platelet simulation. Based on similar trials with stationary objects that were 
roughly perturbed versions of idealized objects in [24], we surmise that the insensitivity to 
the shape parameter in the RBF-IB method stems from the object being deformed in a non- 
smooth fashion by both IB methods. If one were to use our RBF representation within an IB 
method that guaranteed smooth deformations of the immersed elastic objects, we recommend 
choosing as small a shape parameter as the conditioning of the RBF interpolation matrix 
allows. This was shown to give spectral geometric accuracy in [24]. 



6.3 Energy estimates 
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Fig. 5. Energy in the RBF-IB method. The figure on the left shows the energy change per time-step 
computed in the elasticity test as a function of time. The figure on the right shows the energy 
change per time-step computed in the three-platelet simulation as a function of time. For each 
platelet, the number of data sites was set to A'd = 27 and the number of sample sites to A's = 100. 
Both tests were run on a 64 x 128 grid with At = 10~^. 

In this section, we compute energy estimates for the RBF-IB method in two settings: the 
elasticity test and a simulation involving three platelet-like objects, i.e., objects that are 
initially elliptical. The elasticity test is initialized as before and run for 10^ time-steps. The 
platelet simulation is initialized with a parabolic velocity profile and is subject to constant 
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forcing. We halt the simulation when the platelets leave the domain. Both settings use a 
64 X 128 grid and a time-step of At = 10"'^. 

In the elasticity test, one expects the changes in energy to be mainly due to the deformation 
of the stiff elastic object. Eventually, the energy of the system must damp out to zero. The 
platelet simulation, on the other hand, presents a more complex scenario. Depending on 
where the platelets are positioned along the height of the channel, different platelets will 
leave the channel at different times. Also, different platelets will deform to different extents 
depending on where they are in the channel. In this situation, it is unclear what the energy 
profile should look like. Nevertheless, we compute the energy to demonstrate that the energy 
is bounded within the RBF-IB simulation. 

The energy change in a time-step is computed as the sum of the difference in kinetic energy 
of the fluid over the time-step and the potential energy of the elastic object. This can be 
written as 



/ pw"+i ■ n"+^ - / pn" ■ w" + At / F ■ ^ (43) 

J fluid J fluid Jx Ot 



Here, the Lagrangian force F is computed at time level n + 1 at the sample sites. The ^ 
term is computed by applying the evaluation matrix 8^ to the velocities obtained at the data 
sites. This gives us sample site velocities, allowing us to compute dot products with the F 
terms. We compute the discrete analog of Equation (43) for both the elasticity test and the 
platelet simulation. 

The results of this test are shown in Figure 5. The plot on the left shows the change in energy 
of the system for the elasticity test and the plot on the right shows the change in energy of the 
system for the three-platelet simulation. In the elasticity test, the fluid starts off stationary, so 
the initial kinetic energy is zero. However, the elliptical elastic object starts off under tension, 
since its reference configuration is a circle. This means that the initial elastic potential energy 
of the system is high. As the elastic object attempts to minimize its elastic potential energy, 
its deformation drives a change in the kinetic energy of the fluid, causing the kinetic energy 
of the fluid to increase from its initial value of zero to some maximum. However, the elastic 
object soon attains something close to its reference configuration, causing the kinetic energy 
of the fluid to drop sharply. The viscosity of the fluid causes the kinetic energy to eventually 
damp out almost completely, with minor perturbations due to possible deformations of the 
elastic object. The energy of the system continues to decrease as the object becomes more 
and more circular. 

The plot on the right side of Figure 5 shows the change in energy per time-step for the three- 
platelet simulation. The large dip corresponds to two platelets getting highly deformed at 
the edge of the domain, with the increase from that minimum corresponding to the platelets 
leaving the domain. The second dip corresponds to the last platelet leaving the domain. With 
A'd = 27 and A's = 100, the simulation was stable and the energy was bounded. Indeed, the 
RBF-IB method can be stably used to simulate the behavior of over 100 platelets in Poiseuille 
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flow. We present timings for these scenarios furtlier on in tliis section. 



6.4 Platelet modeling comparison 



In this section, we present differences between the RBF-IB method and PL-IB method for 
simulations involving platelets, both in terms of differences in sample site positions and 
differences in computational cost. 



6.4- i Differences between platelets simulated by RBF-IB method and PL-IB method 
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Fig. 6. Normalized differences between the platelets simulated by the RBF-IB method and the 
PL-IB method. The figure on the left shows the normalized maximum differences in centers of mass 
of the platelets just before they leave the channel vs the number of data sites. Similarly, the figure 
on the right shows the normalized maximum differences in sample point positions according to 
Equation (42). The shape parameter was set to e = 3.1. 

Here, we examine the differences between the RBF-IB method and the PL-IB method in the 
context of a small simulation, with three platelet-like objects, i.e., objects that are initially 
elliptical. The difference from the Drift Test, of course, is that the platelets interact with the 
fluid and attempt to retain their original shape. While the simulation of three platelets is 
certainly not a realistic scenario, it provides a simple way to pin down differences between 
the two methods. The platelets were simulated under Poiseuille flow. The simulations were 
stopped when all three platelets departed from the channel {t = 1.2 units). We set the 
maximum velocity of the fluid to be lower in this test than in the platelet simulation test in 
the previous section. This allows us to observe the platelets for a little longer than before. 

Figure 6 shows the differences in the sample site locations of the platelets from the RBF-IB 
and PL-IB methods from these tests. We see that as the number of data sites is increased 
from Nd = 8, the differences in sample site positions that result from RBF-IB simulations and 
the PL-IB simulations are reduced. However, as the number of data sites is further increased, 
these reductions do not continue, demonstrating that the methods are not equivalent to one 
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another in the hmit of A'"d = Ng. However, it is worth noting two things: one, the differences 
between the sample site positions of the two methods are nevertheless quite small, and two, 
the PL-IB method does not represent a gold standard. 
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Fig. 7. Normalized RMS differences between data sites in the RBF-IB method and the corresponding 
sample sites in the PL-IB method as a function of time. The differences using a modified version 
of Equation (41), with the sum being over the data sites instead of the sample sites. The figures 
show this difference on a 32 x 64 grid and a 64 x 128 grid. The figure on the left is for N^ = 25 and 
the figure on the right is for A'^j = 50. The shape parameter was set to e = 3.1. 

While the previous test compares sample sites in the RBF-IB method against IB points 
in the PL-IB method, we thought it useful to compare data sites in the RBF-IB method 
against the IB points in the PL-IB method, thereby giving us a direct comparison of the 
points that are updated using fluid velocities in each method. For this test, both methods 
used A'^^ = 100. The RBF-IB method used A^d = 25, a subset of the 100 sample sites used in 
the PL-IB method. Figure 7 shows the differences between data sites in the RBF-IB method 
and the corresponding sample sites in the PL-IB method. The difference is smaller on the 
finer grid than on the coarse grid. The differences arise due to the changes introduced by our 
algorithm (spreading from sample sites that were evaluated from the RBF representation of 
the boundary and interpolating to data sites). The spikes in the graphs correspond to when 
a platelet exits the domain in a method. 



6.4-2 Timings 

We compared runtimes for the RBF-IB method and the PL-IB method for a range of grid 
sizes, numbers of sample sites N^ and numbers of platelets being simulated. We present 
timings for larger time-step sizes with the RBF-IB method. We also present the timings for 
realistic simulations with both methods. For all tests, the shape parameter was set to e = 3.1. 
Profiling was done with gprof. Matrix-vector multiplications and matrix-matrix multiplica- 
tions were performed with DGEMM (BLAS3). The Intel Fortran Compiler (IFORT) with 
an optimization level of 02 was used. 

Timing comparisons as a function of the number of platelets: 
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Fig. 8. Average time per time-step per grid cell width for 10^ time-steps of each simulation method 
as a function of the number of platelets. In the first row, the figure on the left shows timings on a 
32 X 64 grid and the figure on the right on a 64 x 128 grid. The figure below shows timings on a 
128 X 256 grid. The time-step was set to At = 10~^ for the figures on the top row, and was set to 
At = 10~^ for the figure on the bottom. 

Figure 8 shows timings for three grid sizes for each method as a function of the number of 
platelets (A'p) being simulated. The number of sample sites was fixed at Ng = 100 for both 
methods and the number of data sites for the RBF-IB method was set to A^d = 27 for one 
set of tests and then to TVd = 54 for the next set. For the coarsest grid, profiling shows that 
the cost of platelet operations dominates, while the fluid solver dominates on the finer grids. 

However, while in all cases the cost of platelet operations increases as we increase the number 
of platelets, this cost is much lower for the RBF-IB method due to the RBF representation. 
For example, for A^p = 60, the PL-IB method directly spreads forces from, interpolates to, and 
moves a total of 6000 IB points while the RBF-IB method computes forces at 6000 points, 
and interpolates velocities to (and moves) only 1620 points. If the number of platelets is 
doubled to A^p = 120, the PL-IB method now computes forces at, spreads from, interpolates 
to, and moves 12000 points, whereas the RBF-IB method computes forces at 12000 points, 
but interpolates velocities to and moves only 3240 points. 

As we move to finer grids, profiling shows that the cost of the fluid solver takes the larger 
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share of the total cost per time-step. However, the RBF-IB method appears to require fewer 
iterations of the muhigrid solver to converge to the solution (on an average, one fewer 
iteration per time-step than the PL-IB method) and therefore spends less time in the fluid 
solver, thereby giving us another speedup over the PL-IB method. 

We hypothesize that this may be caused by the RBF-IB method spreading smoother forces 
than the PL-IB method. The cost of the RBF-IB method shows better than linear scaling 
with respect to the number of platelets on all the tested grid sizes for these reasons. The 
difference in scales between the costs on the three grids is due to the increasing cost of the 
fluid solver as the grid becomes finer. 

Timing comparisons as a function of the number of sample sites: 
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Fig. 9. Average time per time-step per grid cell width for 10^ time-steps of each simulation method 
as a function of the number of sample points. In the first row, the figure on the left shows timings 
on a 32 X 64 grid and the figure on the right on a 64 x 128 grid. The figure below shows timings on 



a 128 X 256 grid. The time-step was set to At 
to At = 10^^ for the figure on the bottom. 



10 for the figures on the top row, and was set 



The timings in Figure 9 are presented as a function of the number of sample sites for A'p 



60 



and A'^p = 120 platelets respectively. As before, the difference in scales between the various 
grid sizes arises from the increased cost in the fluid solver as the grid becomes finer. As 
in the previous test, profiling shows that the cost of platelet operations dominates on the 
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coarsest grid, while the cost of the fluid solver dominates on the finer grids. The RBF-IB 
method requires interpolation of velocities to (and position updates of) only the data sites, 
whereas the PL-IB method operates on the sample sites for all platelet operations. This 
means that while the cost of computing forces in the RBF-IB method increases (like in the 
PL-IB method) as the number of sample sites is increased, the cost of updating the data 
sites stays fixed, as does the cost of parametrically interpolating the data site positions. In 
the PL-IB method, both these costs increase with an increasing number of sample sites. 
The iterative fluid solver also converges faster for the RBF-IB method than for the PL-IB 
method, contributing to the sharper rise in cost associated with the PL-IB method on the 
finer grids. This explains the almost-weak scaling with respect to the number of sample sites 
that we see with the RBF-IB method on all grids. 

On the coarsest grid, since fluid solver costs are secondary to the cost of platelet operations, 
we see a gentler increase in cost with the PL-IB method than on the finer grids, although the 
difference in cost between the two methods is nevertheless large. Note that Figure 9 shows 
that both methods have the same increase in cost from 60 to 120 platelets as the timings from 
Figure 8 imply. Furthermore, when N^ = 54 (data not presented), we see a similar increase 
in cost as seen in Figure 8. It is worthwhile to point out that since the spatial configuration 
of the sample sites in this test is different from that in the previous test, the timings need 
not (and do not) match up exactly between the two tests for a given total number of sample 
sites in the domain. 

Timings with a large time-step for the RBF-IB method: 
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Fig. 10. Average time per time-step per grid cell width for 15,583 time-steps of the RBF-IB method 
at Ai = 6.5 X 10"'* as a function of the number of platelets. The figure on the left shows timings 
for a 32 X 64 grid and the figure on the right shows timings for a 64 x 128 grid. The number of 
sample sites was set to N^ = 100. For convenience, we have reproduced the timings from Figure 8 
for those grid sizes. 

As mentioned previously, the RBF-IB method is able to take much larger time-steps than 
the PL-IB method. In Figure 10, we present timings with At = 6.5 x 10~^ as a function of the 
number of platelets being simulated. To replicate the total simulation times in the previous 
tests (10 sec), this simulation was run for 15,583 time-steps. We see that each simulation step 
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using the larger time-step takes roughly twice as long, though the scaling with the number 
of platelets is much the same as before. Profiling shows that this is because the multigrid 
solver appears to take more iterations to converge for the larger time-step, i.e., more work is 
done per time-step. This is to be expected, since the sample sites will have moved further in 
the larger time-step. The advantage is that we need to take far fewer time-steps (less than 
16% the original number) for the same simulation time. Clearly, this results in a speedup. 



6. 5 2D platelet aggregation simulation 





Fig. 11. Platelet aggregation simulations achieved with both the traditional and RBF-IB methods. 
The figure on top shows a zoomed-in snapshot of a platelet aggregation simulation achieved with 
the PL-IB method with a time-step of At = 2 x 10~^, the largest stable time-step allowed by the 
PL-IB method for this grid size; the figure on the bottom shows a zoomed-in snapshot of a platelet 
aggregation simulation achieved with the RBF-IB method with a time-step of Ai = 6.5 x 10~^, 
the largest stable time-step allowed by the RBF-IB method for this grid size. Both snapshots were 
taken at simulation time t = 2.6. Both simulations were run on a 64 x 128 grid with the same initial 
configurations. The arrows in both figures show the magnitude and direction of the velocity field. 

We now present the results of 2D platelet aggregation simulations accomplished using the 
RBF-IB method and the traditional IB method. Unlike the other tests in this section, this 
requires inter-platelet and platelet-wall bonding. These bonds are modeled as outlined in 
Section 4, while the process of platelet activation is modeled as outlined in Section 1.2. We 
start our simulation with platelets placed near the wall to easily demonstrate aggregation. 
The unperturbed fluid velocity has a parabolic profile. The results are shown in Figure 11. 
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It is easy to see from the figure that there are no apparent visual differences between the 
aggregates formed by the traditional IB simulation and the RBF-IB simulation, despite the 
changes introduced by the RBF-IB algorithm. 



6. 6 3D platelet simulation 




Fig. 12. The figure shows the pathUnes of marker particles as they flow past the two bound platelets 
in the 3D simulation. For clarity, we show the pathlines on a single slice of the domain. The platelet 
on the bottom is bound to the wall. The simulation was run on a 64 x 64 x 64 grid with a time-step 
of At = 10~6. 

We also performed a 3D platelet simulation using the RBF-IB method. While we do not 
demonstrate aggregation as in the 2D case, we still demonstrate inter-platelet bonds and 
platelet-wall bonds. We place a single platelet near the wall so that it binds to the wall at 
the beginning of the simulation. We then allow position another platelet at the inlet end at 
an appropriate height to be carried close to the first platelet by the flow. The platelets bond 
with each other and the flow diverts around each platelet, as shown in Figure 12. As in the 
2D simulations, we allow several bonds between platelets. 



7 Discussion 



The IB method, as a numerical methodology for applications involving fluid structure in- 
teractions, naturally lends itself to our problem of interest: simulating platelet aggregation 
during blood clotting. In this application, platelets are modeled as immersed elastic struc- 
tures whose shapes change dynamically in response to blood flow and chemistry. In previous 
work [24], we discussed several geometric representations for platelets and compared them 
to the representation used within the traditional IB method. We concluded that an RBF 
geometric model for platelets would prove advantageous in several ways. In this work, we 
have explored the ramifications of using the RBF geometric model within the IB method. 
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and compared the behavior of this new RBF-IB method against that of the traditional IB 
method in the context of platelet aggregation. We have discussed the issue of selecting an ap- 
propriate shape parameter for the RBF-IB method. We have presented comparisons between 
the two methods in terms of differences in positions between Lagrangian markers on platelets 
simulated by both methods (given identical initial conditions). We have also compared the 
computational costs incurred by both methods in the context of platelet simulations. We 
have compared the volume conservation properties of both methods and also the time-step 
restrictions on both. We have also remarked on the energy properties of our method and 
commented on its stability. 

We advocate the adoption of the new RBF-IB method for the following reasons: 

• it facilitates the evaluation of constitutive models of elasticity in a straightforward way; 

• it reduces the volume loss exhibited by the IB method for problems involving very high 
stiffnesses for the Lagrangian objects when the appropriate number of data sites is chosen; 

• it allows for larger time-step sizes than those allowed in the traditional IB method for a 
given grid size; 

• it is more computationally efficient than the traditional IB method, both due to the uti- 
lization of a small number of data sites and due to smoother forces being spread into the 
fluid resulting in a faster convergence from the fluid solver. 

We note that while the RBF-IB method indeed exhibits lower volume loss (for the appro- 
priate number of data sites) than the traditional IB method, the volume loss has not been 
completely eliminated in applications with very high stiffnesses. We will study ways to alle- 
viate (or eliminate) this problem using the RBF paradigm. Further, while our method has 
been observed to be stable, we note that we do not meet the sufficiency condition outlined 
in [20]. We hope to explore techniques to meet this sufficiency condition without losing the 
advantages afforded to us by the RBF-IB method. Finally, an issue with the RBF-IB method 
is that it is dependent on the parametrization of the immersed elastic objects. For objects 
that are not easily parameterized in terms of circles and ellipses, the use of the RBF model 
as presented in our work (wherein the RBFs are restricted to the circle/sphere) may not be 
ideal. In the future, we thus hope to explore the use of RBFs in a meshfree variational form 
within the IB method so as to be able to easily evaluate constitutive models on arbitrary 
shapes. 
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